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CALCULATION  OF  AIRBLAST  FROM  UNDERGROUND 
AMMUNITION  STORAGE  MAGAZINES 

I . INTRODUCTION 

The  storage  of  munitions  in  underground  facilities 
provides  a potential  hazard  to  surrounding  buildings  or 
populations.  Current  placement  of  such  facilities  relative  to 
above-ground  structures  is  based  on  long  standing  curves  and 
scaling  relations  for  the  peak  airblast  overpressure  as  a 
function  of  weight  of  explosives  and  distance  and  angle 
relative  to  the  tunnel  opening.  The  suggested  scaling  relation 
for  a 50  mbar  safety  criterion  is 


o 

where  d is  the  distance  to  50  mbar  (m) 

D is  the  main  tunnel  diameter  (m) 

PQ  the  gas  pressure  in  the  main  passageway  (bars) 
and  F is  a tabulated  directional  factor  . 

This  relationship  is  based  on  a variety  of  explosive 
weight-to-volume  ratios.  However,  recent  storage  designs  call 
for  larger  amounts  of  explosives  and  lower  explosive  storage 
densities  than  were  considered  in  developing  the  relationship. 

It  has  not  been  demonstrated  that  the  same  relation  holds  for 
these  new  conditions.  In  addition,  the  high  cost  of  real 
estate  suggests  that  small  changes  in  the  safe  distances  will 
have  a significant  financial  impact  on  the  storage  of  munitions. 

A coordinated  program  of  experiments  and  calculations 
is  being  conducted  by  the  Norwegian  Defence  Construction 
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Service  (UDCS)  to  reduce  the  uncertainties  in  determination  of 
safe  distances  and  to  define  the  angular  distribution  of 
overpressure.  As  a part  of  this  program,  S-CUBED  provided 
hydrodynamic  flow  calculations  of  the  detonation,  internal 
propagation  and  external  expansion  of  blast  waves  in  a two- 
dimensional  model  designed  to  simulate  the  NDCS  experiment, 
"Test  1".  The  calculation  was  done  in  four  consecutive  phases, 
which  are  described  in  the  next  four  sections.  Codes  used  were 
HULL,  a two  or  three  dimensional,  state-of-the-art,  finite 
difference  Euierian  hydrodynamic  flow  program;  SAP  (Spherical 
Air  PUFF),  a one-dimensional,  spher ically-symmetric  version  of 
HULL;  and  NLAWS,  an  acoustic  wave  propagation  code.  Using 
multiple  phases  in  this  way,  efficiency  of  computer  time  is 
maximized  while  essential  features  of  the  result  are  retained. 
This  summary  report  provides  a description  of  the  calculation 
and  presents  the  results  that  were  obtained. 
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II.  PHASE  1:  DETONATION  AND  INTERIOR  PROPAGATION 


The  first  phase  of  the  calculation  included  detonation 
of  an  explosive  in  an  enclosed  region  which  was  designed  to  simu- 
late the  interior  tunnel  complex  of  the  NDCS  test  model.  Propagation 
of  blast  waves  in  the  enclosed  chambers  and  exit  of  blast  pressure 
at  the  tunnel  mouth  were  monitored  from  time  zero  (initiation  of 
the  detonation)  to  4.24  msec  after  detonation. 

The  calculational  configuration  used  for  this  phase  is 
illustrated  in  Figure  1.  The  calculational  mesh  is  a rectangular 
grid  in  two-dimensions.  The  HULL  code  was  used  for  this  phase,  as 
well  as  for  Phase  2.  In  order  to  simulate  the  three-dimensional 
test  chamber  in  two  dimensions,  it  was  necessary  to  change  the 
orientations  of  the  detonation  and  expansion  chambers.  In  the 
test  configuration,  the  entranceway  and  expansion  chamber  are 
oriented  so  that  their  long  axes  are  horizontal  and  perpendicular 
to  each  other.  The  explosive  chamber  axis  is  vertical.  In  the 
calculational  configuration,  the  right  angle  orientations  are 
retained,  but  all  chambers  lie  in  the  same  plane. 

Also  because  the  coordinate  system  is  rectangular  while 
some  of  the  test  chambers  have  circular  cross  sections,  it  was 
necessary  to  adjust  the  dimensions  slightly  in  order  to  maintain 
the  appropriate  cross-sectional  areas  and  chamber  volumes.  If  the 
unit  depth  (into  the  paper  in  Figure  1)  is  assumed  to  be  5 cm, 
then  all  of  the  chamber  cross  sections  and  volumes  are  the  same  as 
for  the  test  configuration.  The  total  internal  volume  of  the  complex 
is  17,212  cm3.  Dimensions  of  the  undetonated  explosive  were  also 
adjusted  so  that,  using  the  same  5-cm  depth  and  a density  of 
1.66  gms/cm3 , the  mass  of  the  charge  was  200  gms , corresponding 
to  that  used  in  the  test. 

Grid  size  for  this  calculation  is  approximately  0.5  cm 
vertically  by  0.8  cm  horizontally,  although  adjustments  were  made 
so  that  cell  boundaries  would  coincide  with  the  previously  deter- 
mined material  boundaries.  There  are  70x376,  or  approximately 
26,000,  calculational  cells  in  the  grid. 
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Storage  Bunker  Calculation  (not  to 


The  calculation  was  begun  at  a time  set  arbitrarily 
to  10  ysec.  At  this  time,  the  calculational  zone  at  the 
right-hand  edge  of  the  explosive  was  considered  to  be  deton- 
ated. The  detonation  wave  progressed  through  the  explosive 
from  right  to  left  as  the  calculation  proceeded. 

Two  types  of  graphical  output  are  routinely  provided 
by  HULL,  and  selected  examples  are  included  in  this  report. 

The  first  is  "contour"  plots,  in  which  isograms  of  any  of  the 
hydrodynamic  variables  are  shown  throughout  the  region  of 
interest.  Figures  2 through  9 are  contour  plots  of  pressure 
and  energy  at  four  different  times  during  the  interior  calcu- 
lation. In  Figures  2 and  3,  at  150  ysec,  it  can  be  seen  that 
the  explosive  has  detonated  but  the  energy  has  not  yet  started 
to  escape  from  the  explosive  chamber.  Pressure  and  energy 
values  throughout  the  chamber  are  very  high. 

In  Figures  4 and  5,  at  340  ysec,  a shock  wave  has 
traveled  through  the  narrow  passageway,  reflected  from  the 
left-hand  side  of  the  vertical  chamber,  and  is  moving  down- 
ward. Pressures  are  very  high  at  the  left-hand  side  of  the 
vertical  chamber  where  the  flow  has  stagnated,  and  reflected 
shocks  from  the  corners  are  visible  in  the  explosive  chamber. 

The  next  two  figures,  6 and  7,  illustrate  the  situa- 
tion of  500  ysec.  At  this  time,  energy  is  moving  down  the 
long  entrance  tunnel.  The  interesting  point  about  the  flow 
in  this  region  is  that  the  energy  contours  are  not  perpendi- 
cular to  the  tunnel  walls.  Because  of  the  right-angle  bend, 
energies  are  higher  near  the  bottom  of  the  tunnel. 

The  last  set  of  contour  plots  was  made  at  900  ysec, 
and  is  shown  in  Figures  8 and  9.  In  this  set,  the  shock  has 
reached  the  tunnel  entrance  and  is  expanding  as  it  moves 
into  the  exterior  space. 
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Figure  2.  Pressure  Contours  in  Explosive  Chamber 
at  150  usee  after  Detonation 
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Figure  3.  Energy  Contours  in  Explosive  Chamber  at 
150  usee  after  Detonation 
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Figure  4.  Pressure  Contours  in  Vertical  Chamber  and 
Passageway  at  340  ysec  after  Detonation 
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Figure  5.  Energy  Contours  in  Vertical  Chamber  and 
Passageway  at  340  ysec  after  Detonation 
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Figure  6.  Pressure  Contours  in  Vertical  Chamber  and 

Entrance  Tunnel  at  500  usee  after  Detonation 
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Figure  7.  Energy  Contours  in  Vertical  Chamber  and  Entrance 
Tunnel  at  500  usee  after  Detonation 
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Figure  8 . Pressure  Contours  in  Tunnel  Entrance  Region 
at  900  nsec  after  Detonation 
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The  second  type  of  output  is  probably  of  more  imme- 
diate interest  in  this  application  because  it  is  directly 
comparable  to  experimental  results.  "Stations"  are  predeter- 
mined locations  in  the  calculational  grid  at  which  hydro- 
dynamic  data  are  monitored  as  functions  of  time.  For  the 
interior  phase  calculation,  stations  were  placed  at  measuring 
points  corresponding  to  those  in  the  NDCS  test.  These  are 
the  points  labeled  1 through  5 in  Figure  1.  Additional 
stations  were  positioned  at  the  tunnel  entrance  and  at  other 
points  throughout  the  calculational  grid.  Stations  28,  29 
and  30  are  of  particular  interest  because  they  are  the  ones 
used  to  drive  the  second  phase  of  the  calculation.  Their 
locations  are  also  shown  in  Figure  1. 

Figures  10  through  14  are  overpressure  versus  time 
records  at  Stations  1 through  5.  The  first  three,  which  are 
in  the  explosive  and  vertical  chambers  and  the  connecting 
tunnel  between  them,  show  a large  initial  spike  followed  by 
many  smaller  peaxs.  The  latter  are  due  to  reflections  from 
the  chamber  walls.  Stations  4 and  5,  for  which  the  over- 
pressure records  are  shown  in  Figures  13  and  14,  are  located 
in  the  entrance  tunnel.  The  records  show  many  sharp  spikes, 
which  are  a result  of  the  mixture  of  air  and  detonation 
products  flowing  in  this  region.  The  maximum  overpressure 
in  the  entrance  tunnel  is  about  80  bars.  At  Measuring  Point 
5,  the  density  versus  time  record  is  given  in  Figure  15. 
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Figure  11.  Overpressure  versus  Time  Record  at 
Measuring  Point  2 
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Figure  13.  Overpressure  versus  Time  Record  at 
Measuring  Point  4 
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Parameters  at  the  tunnel  entrance  (Stations  28,  29 
and  30)  are  shown  in  Figures  16  through  21.  For  these  stations, 
the  horizontal  flow  velocities  versus  time  are  given  as  well 
as  the  overpressures.  Note  that  these  velocities  are  negative 
because  flow  through  the  portal  is  in  the  "-x"  direction. 
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Figure  17.  Horizontal  Velocity  versus  Time  Record  at 
Tunnel  Entrance,  Near  the  Floor 
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Figure  18.  Overpressure  versus  Time  Record  at  Tunnel 
Entrance,  Near  the  Center 
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Figure  19.  Horizontal  Velocity  versus  Time  Record  at 
Tunnel  Entrance,  Near  the  Center 
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Figure  20.  Overpressure  versus  Time  Record  at  Tunnel 
Entrance,  Near  the  Roof 
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Figure  21.  Horizontal  Velocity  versus  Time  Record  at 
Tunnel  Entrance,  Near  the  Roof 
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III.  PHASE  2:  EXPANSION  INTO  EXTERIOR  REGION 


In  order  to  model  the  exterior  region,  a computational 
mesh  was  set  up  using  cylindrical  coordinates.  The  mesh  con- 
sists of  182x176  zones,  for  a total  of  31,304.  The  overall 
grid  dimensions  are  5.96  m by  5.99  m.  A 100x100  cell  subgrid, 
in  which  the  size  of  each  cell  is  approximately  0.49x0.8  cm, 
was  defined  on  the  cylindrical  axis  at  the  tunnel  opening. 

Beyond  the  subgrid,  cell  size  is  expanded  by  anout  5 percent 
per  cell  in  each  direction.  Stations  were  located  at  intervals 
along  radial  lines  from  the  tunnel  exit. 

The  configuration  is  illustrated  in  Figure  22-.  For 
this  phase,  a plane  through  the  axis  of  symmetry,  which  can 
be  thought  of  as  a flat,  perfectly  reflecting  surface,  becomes 
the  ground  plane.  The  tunnel  opening  and  cliff  face  become 
cylindrical  sections  when  converted  from  Cartesian  to  cylin- 
drical coordinates.  The  radius  of  the  opening  was  adjusted  so 
that  the  cross-sectional  area  of  the  naif-disk  opening  in 
cylindrical  coordinates  is  equal  to  the  rectangular  5x6  cm 
tunnel  opening  of  Phase  1. 

The  cliff  face  was  modeled  by  placing  a row  of  "island" 
cells  along  the  bottom  of  the  mesh.  An  island  is  a non-compres- 
sible,  reflective  cell  of  the  same  type  as  were  used  in  Phase  1 
for  the  tunnel  walls.  The  configuration,  it  should  be  noted, 
is  rotated  by  90  degrees  from  that  of  Phase  1 , so  that  direc- 
tion of  flow  is  upward  in  Figure  22. 

The  calculation  was  initiated  at  0.7  msec  by  setting 
the  hydrodynamic  parameters  in  the  cells  at  the  tunnel  opening 
to  those  read  from  the  station  records  of  Phase  1 . To  perform 
the  transformation,  the  "-x"  velocities  of  the  Phase  1 station 
records  became  the  "+y"  velocities  of  the  Phase  2 driver  cells. 


28 


and  +y ' s became  +x's.  Interpolation  was  performed  to  drive 
the  cells  at  locations  between  stations  at  appropriate  values. 
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Figure  22.  Calculational  Configuration  for  Phase  2 


The  second  phase  calculation  was  run  to  3.1  msec  after 
detonation  of  the  explosive.  By  this  time,  the  shock  wave 
emerging  from  the  tunnel  opening  had  traveled  a distance  of 
1.75  m (5.74  ft)  along  the  axis  and  had  decayed  to  an  over- 
pressure of  approximately  1.0  bar  (14.5  psi)  . Figures  23 
through  28  illustrate  the  pressure,  density,  and  energy  at 
two  different  times  (1.5  and  2.3  iyisec)  during  the  calculation. 
An  interesting  feature  shown  by  the  pressure  and  density  plots 
is  the  development  of  a Mach-stem-like  structure  at  the  inter- 
action with  the  cliff  face.  The  structure,  usually  observed 
only  with  height-of-burst  detonations , is  an  indication  of 
flow  toward  the  surface  (downward  in  the  picture)  as  well  as 
parallel  to  the  surface. 
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Figure  23 . Pressure  Contours  in  Exterior  Region  at 
1.5  msec  after  Detonation 
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Figure  24 . Density  Contours  in  Exterior  Region  at 
1.5  msec  after  Detonation 
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Figure  26.  Pressure  Contours  in  Exterior  Region  at 
2.3  msec  after  Detonation 
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Figure  27.  Density  Contours  in  Exterior  Region  at 
2.3  msec  after  Detonation 
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Figure  28.  Energy  Contours  in  Exterior  Region  at 
2.3  msec  after  Detonation 
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IV.  PHASE  3:  ONE-DIMENSIONAL,  CALCULATIONS 


Phase  3 consists  of  a series  of  one-dimensional 
calculations  along  lines  radiating  from  the  tunnel  entrance. 

The  Phase  3 calculations  were  done  using  SAP  (Spherical  Air 
Puff),  S-CUBED 1 s version  of  the  one-dimensional  Eulerian, 
spherical- symmetry  hydrocode  originally  developed  at  the  Air 
Force  Weapons  Laboratory  (AFWL) . 

At  3.1  msec,  as  previously  mentioned,  the  Phase  2 
two-dimensional  calculation  was  terminated.  It  was  replaced  by 
a series  of  one-dimensional  calculations  along  the  lines 
A through  G shown  in  Figure  29.  The  initial  conditions  for 
each  Phase  3 calculation  were  obtained  by  "cutting"  through  the 
Phase  2 calculational  mesh  at  3.1  msec.  Energy,  density  and 
velocity  values  were  read  from  zones  through  which  the  lines 
passed.  Interpolation  across  four  adjacent  zones  was  performed 
in  order  to  obtain  values  along  the  lines.  Also,  the  component 
x and  y velocities  were  combined  to  determine  the  appropriate 
velocity  in  the  direction  of  the  line  cut.  Phase  2 contour 
plots  at  3.1  msec,  which  were  used  to  provide  the  initial 
conditions,  are  given  in  Figures  30  through  33.  The  lines 
along  which  the  cuts  were  taken  (approximately  every  15  deg) 
are  shown  in  Figure  30. 

As  each  of  the  SAP  calculations  proceeded  in  time,  the 
waveform  defined  by  the  cut  was  allowed  to  propagate  outward. 
The  innermost  cell  was  continually  reset  to  conditions  from  the 
lat6-time  record  tor  Station  29  (Phase  1).  The  one-dimensional 
calculations  were  continued  until  a shock  front  peak  over- 
pressure of  less  than  250  mbar , or  3.62  psi,  was  reached.  This 
occurred  at  times  between  12.0  and  12.5  msec.  After  4.2  msec, 
when  the  Station  29  record  ended,  constant  values  of  the 
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parameters  equal  to  those  at  the  last  time  for  Station  29 
(4.2  msec)  were  maintained  at  the  inner  edge  of  the  SAP  grid. 


Figure  29.  Calculational  Conf icruration  Showina 
Lines  Along  Which  One-Dimensional 
Calculations  Were  Performed 
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Figure  30.  Pressure  Contours  at  3.1  msec  Used  for 
Initiation  of  Phase  3 
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Figure  31.  Energy  Contours  at  3.1  msec  Used  for 
Initiation  of  Phase  3 
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Figure  32.  Density  Contours  at  3.1  msec  Used  for 
Initiation  of  Phase  3 
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Figure  33. 


Flow  Velocity  Vector  Plot  at  3.1  msec  Used  for 
Initiation  of  Phase  3 
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Figures  34  through  40  are  overpressure  versus  range 
plots,  one  for  each  of  the  seven  SAP  calculations,  at  the  final 
time  for  that  calculation.  These  range  plots  provided  the 
waveforms  from  which  inputs  for  NLAWS , the  acoustic  wave 
propagation  code  used  in  Phase  4,  were  defined. 
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Figure  34.  Overpressure  versus  Range  Plot  Along  90° 
Line  at  12.15  msec 
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Figure  35-  Overpressure  versus  Range  Plot  Along  75°  Line 
at  12. QO  msec 
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Figure  36.  Overpressure  versus  Range  Plot  Along  60°  Line 
at  12.00  msec 
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Figure  37.  Overpressure  versus  Range  Plot  Along  45°  Line 
at  12.50  msec 
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Figure  38.  Overpressure  versus  Range  Plot  Along  30°  Line 
at  12.50  msec 
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Figure  39.  Overpressure  versus  Range  Plot  Along  15°  Line 
at  12.50  msec 


48 


OVERPRESSURE  DYNES/CM**2 


S-CUBED 


H 

L_i 


NORWAY  AMMO  BUNKER-EXTERIOR  LINE  F 


Figure  40.  Overpressure  versus  Range  Plot  Along  0°  Line 
at  12.50  msec 
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V.  PHASE  4:  NONLINEAR  ACOUSTIC  WAVE  SOLUTION  (NLAWS ) 


NLAWS  is  a code  which  treats  the  airblast  waveform 
as  an  acoustic  wave.  The  code  assumes  an  ideal  triangular 
waveform  input  with  an  infinitely  sharp  rise  to  a positive 
overpressure  and  then  a monotonic  decay  to  zero.  The  triangu- 
lar wave  is  then  propagated  acoustically  until  it  has  attenuated 
to  the  desired  level.  NLAWS  was  written  at  the  Air  Force 
Weapons  Laboratory  and  is  used  to  find  the  peak  overpressure 
and  velocity  at  ranges  below  approximately  280  rabar  (4  psi) . 

The  NLAWS  program  was  run  for  each  of  the  one- 
dimensional calculations,  beginning  at  the  last  time  from  the 
SAP  calculations  of  Phase  3.  The  objective  was  to  determine 
a range  at  which  the  peak  overpressure  would  have  attenuated 
to  50  mbar.  Because  the  input  waveform  for  NLAWS  must  be 
triangular,  waveforms  from  the  SAP  one-dimensional  calculations 
had  to  be  modified.  This  was  done  in  the  following  manner. 
First,  a vertical  line  was  drawn  through  the  mid-point  of  the 
first  rise  on  the  SAP  waveform.  In  Figure  41,  this  line  is 
designated  "LINE  1".  Second,  a representative  slope  was 
chosen  on  the  falling  portion  on  the  wave,  and  a line  was 
drawn  tangent  to  the  sloping  portion  (illustrated  by  "LINE  2" 
in  Figure  41) . These  two  lines  defined  a peak  overpressure 
and  a positive  phase  duration  for  input  to  NLAWS.  Because 
there  was  some  uncertainty  about  the  choice  of  the  tangent 
line,  a second  line  (LINE  3)  was  chosen  and  a second  NLAWS 
calculation  was  performed.  Choosing  the  input  waveform  in 
this  manner  eliminates  numerical  shock-front  smearing  and 
any  numerical  overshoot  from  the  hydrocode  results.  The  two 
different  NLAWS  calculations  provide  an  estimate  of  the 
accuracy  of  the  results. 


50 


OVERPRESSURE  DYNES/CM**2 


Figure  41.  Overpressure  vs.  Range  Plot  Showing  Method  of 
Determining  Triangular  Wave  Input  for  NLAWS 
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The  procedure  described  previously  was  applied  to 
each  of  the  seven  waveforms  from  the  SAP  calculations  given 
in  the  last  section.  In  cases  such  as  that  for  the  45  deg 
line  (Figure  37) , in  which  waveform  consisted  of  two  distinct 
pulses,  only  the  first  (and  largest)  was  used.  Under  acoustic 
conditions,  all  parts  of  a wave  propagate  at  the  same  speed, 
thus  there  is  no  overtaking  of  tne  front-running  portions  by 
those  following. 
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VI.  CONCLUSIONS 


The  peak-overpressure  versus  range  listings  generated 
by  NLAWS  as  described  in  the  previous  section  were  interpolated 
as  necessary  to  determine  the  range  at  which  this  peak  falls 
below  50-mbars.  Because  two  sets  of  input  conditions  were 
defined,  a spread  in  range  was  obtained  at  each  angle.  The 
results  are  shown  in  Figure  42  (outer  curve).  As  can  be  seen, 
for  positions  directly  in  front  of  the  tunnel  entrance  (0  deg), 
a distance  of  17  meters  is  required  to  assure  that  the  peak 
overpressure  experienced  will  be  less  than  50  mbars.  For 
otf— axis  positions,  16  meters  is  sufficient  at  15  and  30  deg 
and  14  meters  at  45  deg.  At  angles  of  60  deg  or  greater,  only 
11  meters  is  required.  For  the  full-scale  situation  of 
interest,  in  which  200,000  kgm  of  explosive  are  detonated,  all 
distances  would  be  scaled  up  by  a factor  of  100. 

The  inner  curve  in  Figure  42  defines  the  200-mbar  peak 
overpressure  distances.  This  curve  was  determined  in  exactly 
the  same  way  as  the  50-mbar  curve,  except  that  in  this  case  the 
results  were  taken  from  the  early  part  of  the  NLAWS  run  or,  at 
several  of  the  angles,  directly  from  the  SAP  results. 
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Figure  42.  Ranges  for  50  and  200  mbar  Peak  Overpressures 
as  Functions  of  Angle 
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